home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / robust / como.ex next >
Text File  |  1999-09-16  |  10KB  |  465 lines
  1. // Examples of Standard plants (mostly from Kwakernak)
  2. // Define the polynomial s by typing "s=poly(0,'s')"
  3. // Most examples are in transfer representation but computations
  4. // are done in state space form
  5.  
  6. //2.4.1 (Kwakernaak)
  7.  
  8. H=[1 (s-1)/((s-2)*(s-3));
  9.    1 (s-1)/((s-2)*(s-3))];
  10. r=[1 1];
  11.  
  12. //Transform to get d12 full rank
  13.  
  14. H(:,2)=H(:,2)*(s+1);
  15.  
  16. [sk,mu]=h_inf(H,r,0,1,20);
  17.  
  18. //transform back
  19.  
  20. sk=sk*(s+1);
  21.  
  22. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  23.  
  24. //3.9.1 
  25. c=0.1;r=0;
  26.  
  27. H=[(1+sqrt(2)*s+s**2)/(s*s) 1/(s*s);
  28. 0 c+c*r*s;
  29. (1+sqrt(2)*s+s**2)/(s*s) 1/(s*s)];
  30.  
  31. r=[1 1];
  32.  
  33. //[sk,mu]=h_inf(H,r,0,1,20);
  34.  
  35. //compare (185) in como notes
  36.  
  37. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  38.  
  39. c=0.1;r=1;
  40.  
  41. H=[(1+sqrt(2)*s+s**2)/(s*s) 1/(s*s);...
  42.    0 c+c*r*s;..
  43.    (1+sqrt(2)*s+s**2)/(s*s) 1/(s*s)];
  44.  
  45. r=[1 1];
  46.  
  47. // divide 2nd column of h by (s+2) for properness
  48. H(:,2)=H(:,2)/(s+2);
  49.  
  50. //[sk,mu]=h_inf(H,r,0,1,30);
  51.  
  52. // Remove parasitic root introduced above :
  53. //vv=roots(sk(2))   //vv(2) = -2 approximately 
  54.  
  55. //sk=sk/real((s-vv(2)))
  56.  
  57. //compare (186) como notes
  58.  
  59. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  60.  
  61. //3.12.1
  62.  
  63. c=0.1;
  64. H=[(s+1)/s 1/s;
  65.   0  c;
  66.  (s+1)/s 1/s];
  67.  
  68. r=[1 1];
  69.  
  70. //[sk,mu]=h_inf(H,r,0,1,20);
  71.  
  72. //sk=clean(sk,0,0.00001);
  73.  
  74. //compare (228)
  75.  
  76. [sk,mu]=h_inf(H,r,0,10,20);
  77.  
  78. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  79.  
  80. //4.12.1
  81.  
  82. H=[(s*s+sqrt(2)*s+1)/(s*s) 1/(s*s);
  83.  0 1;
  84. (s*s+sqrt(2)*s+1)/(s*s) 1/(s*s)];
  85.  
  86. r=[1 1];
  87.  
  88. //[sk,mu]=h_inf(H,r,0,10,20);
  89.  
  90. //compare (422) which is misprinted !
  91.  
  92. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  93.  
  94. //1990 IFAC (Kwakernaak)
  95.  
  96. H=[-(1+s/4)/(1+2*s) (s+1)/(s-2) 0 1/(s-2) 0;
  97.    0 0 0 1 0;
  98.    0 0 0 0 1;
  99. (1+s/4)/(1+2*s) 0 0 0 0;
  100.  0 (s+1)/(s-2) 1 1/(s-2) 0];
  101.  
  102. r=[2 2];
  103.  
  104. //[Sk,mu]=h_inf(H,r,0,1,30);
  105.  
  106. //compare (33) in Matlab report 
  107.  
  108. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  109.  
  110. //Mc-Farlane Glover for 1/s^2
  111.  
  112. G=[(1+s*sqrt(2)+s*s)/(s*s) 1/(s*s);
  113.    0 1;
  114.    (1+s*sqrT(2)+s*s)/(s*s) 1/(s*s)]
  115.  
  116. r=[1,1]
  117.  
  118. //[sk,mu]=h_inf(G,r,0,1,20)
  119.  
  120. //compare (25) in Matlab report
  121.  
  122. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  123.  
  124. //Mc Farlane Glover for 1/s^3
  125.  
  126.  
  127.  
  128. G=[(1+2*s+2*s*s+s**3)/(s**3) 1/(s**3);
  129.   0 1;
  130.   (1+2*s+2*s*s+s**3)/(s**3) 1/(s**3)]
  131.  
  132. r=[1 1];
  133. //[sk,mu]=h_inf(G,r,0,1,20)
  134.  
  135. //compare (27) in Matlab report
  136.  
  137. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  138.  
  139. //SISO Model matching
  140.  
  141. G=[6*(s-1)/((s-2)*(s-3)) 1;
  142.    -1 0];
  143. r=[1 1];
  144.  
  145. //multiply first row by all-pass transfer to have stabilizability
  146. //and detectability
  147.  
  148. k=((s-2)*(s-3))/((s+2)*(s+3))
  149.  
  150. G(1,:)=k*G(1,:)
  151.  
  152. [sk,mu]=h_inf(G,r,0,5,20)
  153.  
  154. //compare (42) in Matlab report
  155.  
  156. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  157.  
  158. //Mimo Model matching
  159.  
  160. G=[0.5/(s-1) 0 1 0;
  161.    1/(s*s-s+1) 2/(s-1) 0 1;
  162.     -1 0 0 0;
  163.      0 -1 0 0];
  164.  
  165. r=[2 2]
  166.  
  167. p=(s*s-s+1)*(s-1);
  168.  
  169. q=p/horner(p,-s);
  170.  
  171. H=[q 0;
  172.    0 q]
  173.  
  174. //multiplication from the left by inner H
  175.  
  176. G(1:2,:)=H*G(1:2,:);
  177.  
  178. //Stabilizability and detectability assumptions are now verified
  179.  
  180. [sk,mu]=h_inf(G,r,0,1,30);
  181.  
  182. // the parasitic (?) root -83.60 (Kwakernaak) or -84.53 (Francis) is not
  183. // in the poles of sk :
  184. // one finds the roots -0.6920857 + or -0.7932064i and -1.1030557
  185.  
  186. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  187.  
  188. //Mixed sensitivity (ex 11 in matlab notes)
  189.  
  190. G=[(s+1)/s 0 1/s 0;
  191.    0 (s+0.5)/(s-1) 0 1/(s-1);
  192.    0 0 1 0;
  193.    0 0 0 1;
  194. (s+1)/s 0 1/s 0;
  195. 0 (s+0.5)/(s-1) 0 1/(s-1)];
  196.  
  197. r=[2 2]
  198.  
  199. [sk,mu]=h_inf(G,r,0,1,30);
  200.  
  201. sk=clean(sk,1.d-6)
  202.  
  203. //compare (40) in Matlab notes
  204.  
  205. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  206.  
  207. //Laub's ex (CDC 90)
  208.  
  209. A=[0.0904 0.8465;
  210.    0.6888 0.7152];
  211.  
  212. B=[0.7092 0.3017 0.7001;
  213.    0.1814 0.9525 0.1593];
  214.  
  215. C=[0.3088 0.7350;
  216.    0.5735 0.9820;
  217.    0.6644 0.5627];
  218.  
  219. D=[0.7357 0.4156 0.0;
  220.    0.2588 0.1544 1.0;
  221.     0.0    1.0   0.0];
  222.  
  223. r=[1 1];
  224.  
  225. p=syslin('c',a,b,c,d);
  226.  
  227. //[sk,mu]=h_inf(P,r,0,10,20);
  228.  
  229. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  230.  
  231. //Mustafa's ex
  232.  
  233. A=[20 -100;1 0];
  234. B=[1;0];
  235. C=[1,-0.1];
  236. Gss=syslin('c',A,B,C);
  237. [Pss,r]=normlqg(Gss);    //The standard plant
  238.  
  239. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  240.  
  241. //Maciejowski's example
  242.  
  243. At=[-0.001 0 1.1320 0 -1;
  244.  0 -0.0538 -0.1712 0 0.0705;
  245.  0 0 0 1 0;
  246.  0 0.0485 0 -0.8556 -1.013;
  247.  0 -0.2909 0 1.0532 -0.6859];
  248.  
  249. B=[0 0 0;
  250.    -0.120 1 0;
  251.  0  0 0;
  252. 4.4190 0 -1.665;
  253. 1.575 0 -0.0732];
  254.  
  255. C=[1 0 0 0 0;
  256.    0 1 0 0 0;
  257.    0 0 1 0 0];
  258.  
  259. Dt=1.d-5*eye(3,3);
  260. D=0*Dt;
  261.  
  262. Gss=syslin('c',At,b,c,D);
  263.  
  264. //gtf=ss2tf(Gss);gtf=clean(gtf,1.d-10);
  265.  
  266. w1=syslin('c',((s+6)**2)/((s+0.00006)*(s+0.6)));
  267.  
  268. w2=syslin('c',2000*(s+10)*(s+50)/((s+1000)**2));
  269.  
  270. W1ss=tf2ss(w1);
  271. W2ss=tf2ss(w2);
  272.  
  273. W1=[W1ss,0,0;0,W1ss,0;0,0,W1ss];
  274. W2=[W2ss,0,0;0,W2ss,0;0,0,W2ss];
  275.  
  276. A1=W1(2);B1=W1(3);C1=W1(4);D1=W1(5);
  277. A2=W2(2);B2=W2(3);C2=W2(4);D2=W2(5);
  278.  
  279. //Ptf=[W1, W1*Gtf;
  280. //     0*eye(3,3), W2*Gtf;
  281. //     eye(3,3), -Gtf];
  282.  
  283. //Pss=tf2ss(Ptf);
  284.  
  285. r=[3,3];
  286.  
  287. Ap=[At,0*eye(5,6),0*eye(5,6);
  288.    B1*C,A1,0*eye(6,6);
  289.    B2*C,0*eye(6,6),A2];
  290.  
  291. Bp=[0*eye(5,3),B;
  292.    -B1,B1*D;
  293.     0*eye(6,3),B2*D];
  294.  
  295. Cp=[-D1*C,-C1,0*eye(3,6);
  296.     D2*C,0*eye(3,6),C2;
  297.     -C,0*eye(3,6),0*eye(3,6)];
  298.  
  299. Dp=[D1,-D*Dt;
  300.     0*eye(3,3),D2*Dt;
  301.    eye(3,3),-D];
  302.  
  303. P=syslin('c',Ap,Bp,Cp,Dp);
  304.  
  305. [sk,mu]=h_inf(P,r,0,1,30);
  306.  
  307. //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  308.  
  309.  
  310. //Safonov-Chiang
  311.  
  312. A=[-2.2567e-02  -36.617  -18.897  -32.090  3.2509  -7.6257e-1;
  313.   9.2572e-05  -1.8997  9.8312e-01  -7.2562e-04  -1.7080e-01  -4.9652e-03;
  314.     1.2338e-02  11.720  -2.6316  8.7582e-04  -31.604  22.396;
  315.     0           0        1       0            0       0  ;
  316.     0           0        0       0           -30      0  ;
  317.     0           0        0       0            0      -30  ];
  318.  
  319. B=0*ones(6,2);B(5,1)=30;B(6,2)=30;
  320. C=0*ones(2,6);C(1,2)=1;C(2,4)=1;
  321. D=0*ones(2,2);
  322. G=syslin('c',a,b,c);r=syssize(G);
  323. W1=[(s+0.01)/(1+0.01*s)  0;0    (s+0.01)/(0.01*s+1)], 
  324. tau=0.0005;
  325. W2=[1000/(s*s)   0;0    1000/(s*s*(s*tau+1))],
  326.  
  327.  
  328. //
  329.  
  330. s=poly(0,'s');
  331. G=(s+1)/(s*s+1);
  332. w1=1;w3=(s+2);
  333.  
  334. P=[w1, -w1*G;
  335.    0 , w3*G
  336.    1 , -G];
  337.  
  338. r=[1,1];
  339. Pss=tf2ss(P);
  340. P12=Pss(1:2,2);P21=Pss(3,1);P22=Pss(3,2);
  341. sl=p21;
  342. sm=[s*eye-sl(2),sl(3);sl(4),-sl(5)];
  343. detr(sm)       //Poles of G --> Transmission zeros of P21
  344. [sk,mu]=h_inf(P,r,0,1,20);   //Fails
  345.  
  346. Ptmp=P;fil=(1+s*s)/(1+s)/(1+s);Ptmp(:,2)=Ptmp(:,2)*fil;
  347. [sk,mu]=h_inf(Ptmp,r,0,1,10);    //---> sk=0 mu=1 OK
  348.  
  349. // Random example...
  350.  
  351. nx=5;
  352. nu=2;ny=3;
  353. nw=4;nz=3;
  354. A =[1,2,3,0,1;4,5,6,1,0;7,8,9,1,1;0,1,0,1,1;5,4,3,2,1];
  355. B2=[1,2;2,1;1,1;1,0;2,0];  //(nx,nu)
  356. C2=[1,1,0,1,0;0,1,1,2,1;2,3,1,0,1];   //(ny,nx)
  357. r=[ny,nu];
  358.  
  359. //%Q=rand(nx+nu,3)*rand(3,nx+nu);Q=Q*Q';
  360. //%R=rand(nx+ny,2)*rand(2,nx+ny);R=R*R';
  361.  
  362.  
  363. D12=[0,1;2,1;3,0];    //%(nz,nu)
  364. D21=[1,2,1,1;2,1,1,2;0,2,1,1];    //%(ny,nw);
  365. C1=[1,1,0,1,1;2,1,1,0,1;1,0,0,1,1];   //%(nz,nx)
  366. B1=[1,0,1,0;2,3,3,1;3,1,1,0;1,2,1,0;0,1,1,0];  //%(nx,nw)
  367. //%Q#=[C1 D12]'*[C1 D12];
  368. //%R#=[B1;D21]*[B1;D21]';
  369. D11=2*[1,2,1,0;2,2,1,1;0,0,2,1]; //%(nz,nw)
  370. D11=D11;
  371. D22=0*[1,2;2,1;3,1];       //%(ny,nu)
  372. gama=64;
  373. B1=B1/sqrt(gama);
  374. D11=D11/sqrt(gama);
  375. D21=D21/sqrt(gama);
  376. B2=B2*sqrt(gama);
  377. D12=D12*sqrt(gama);
  378. C1=C1/sqrt(gama);
  379. D11=D11/sqrt(gama);
  380. D12=D12/sqrt(gama);
  381. C2=C2*sqrt(gama);
  382. D21=D21*sqrt(gama);
  383. B=[B1,B2];C=[C1;C2];
  384. D=[D11,D12;D21,D22];
  385.  
  386. P=syslin('c',A,B,C,D);
  387.  
  388. [K,mu]=h_inf(P,r,0,10,20);
  389.  
  390. //    Difficult example gamma between 1.6 and 1.8
  391. // h_inf finds 1.671244...  and
  392. // K=[(9.7574+0.2790*s)/(0.023419+s), 0.1346]
  393. // with P1=minreal(P) (state dimension 5 instead of 8 for P)
  394. // get same gamma but a controller of order (4,3)
  395. w=[-1.05033595281261e+00 1 0  -2.72165256144227e+03 0 0 0 0 0 0 0 0 0 0;
  396.    -1.39973553704802e+02 0 0  -2.94468368703638e+06 0 0 0 0 0 0 0 0 0 0;
  397.    0 0  -2.100e+02  -2.250e+04 0 0 0 0 0   3.80748204520352e-01 0 0 0 1;
  398.    0 0 1 0 0 0 0 0 0 0 0 0 0 0;
  399.    0,0,0,0, -4.89284164859002e+03, -7.98004338394794e+06,...
  400.  -4.33839479392625e+09,0 , 4.49639998582548e-02,0,0,0,0,0;
  401.    0 0 0 0 1 0 0 0 0 0 0 0 0 0;
  402.    0 0 0 0 0 1 0 0 0 0 0 0 0 0;
  403.   1.80466256786985e+00 0 0   4.58687395748091e+03 0 0 0,...
  404.   -2.34192037470726e-02 0 0 1  -1 0 0;
  405.    0 0 0 1.25100080458419e+04 0 0 0 0 0 0 0 0 0 0;
  406.    0 0 0 0 0 0 0 0 0   6.000e-01 0 0 0 1.57584459460774e+00;
  407.    6.31631898754447e-01 0 0 1.60540588511832e+03 0 0 0,...
  408.    4.67564402810304e+00 0 0   3.500e-01  -3.500e-01 0 0;
  409.  1.80466256786985e+00 0 0   4.58687395748091e+03 0 0 0 0 0 0 1  -1 0 0;
  410.    0 1 0 0 0 0 0 0 1.79855999433019e+00 0 0 0 1 0];
  411.  
  412. A=W(1:8,1:8);B=W(1:8,9:14);C=W(9:13,1:8);D=W(9:13,9:14);
  413. P=syslin('c',A,B,C,D);
  414. r=[2,1];
  415.  
  416. [K,mu]=h_inf(P,r,0,1,40);
  417.  
  418. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  419.  
  420. //Automotive fuel control
  421. eps=0.001;
  422. s=poly(0,'s');
  423. P=(5.498*s*s+4.007d+2*s-4.444d+5)/((s+eps)*(s**3+93.72*s*s+9.520d+3*s+1.214d+5));
  424. P=P*(18000/5.498);
  425. ro=2;
  426. w2=(s+30)*(s+60)/30/60/10;
  427. w1=20*ro/((s+eps)*(s+20));
  428. G=[w1,-w1*P;
  429.    0 ,w2*P;
  430.    1, -P];
  431. W1ss=tf2ss(w1);
  432. Pss=tf2ss(P);
  433. W2P=W2*P;
  434. //calcul de cw2;
  435. www=clean((s*eye-Pss(2))**(-1)*bp);
  436. mat=[horner(www,-1),horner(www,1),horner(www,10),horner(www,20)];
  437. ww=clean(w2p-1);
  438. f=[horner(ww,-1),horner(ww,1),horner(ww,10),horner(ww,20)];
  439. cw2=f/mat;
  440. W2pss=syslin('c',Pss(2),Pss(3),cw2,1);
  441.  
  442. AGss=[W1ss(2),W1ss(3)*Pss(4);0*eye(4,2),Pss(2)];
  443. BGss=[-W1ss(3),0*ones(2,1);0*ones(4,1),Pss(3)];
  444. CGss=[-W1ss(4),0*ones(1,4);0*ones(1,2),w2pss(4);0*ones(1,2),-pss(4)];
  445. DGss=[0,0;0,1;1,0];
  446. Gss=syslin('c',AGss,BGss,CGss,DGss);
  447. r=[1,1];
  448. [a,b1,b2,c1,c2,d11,d12,d21,d22]=smga(Gss,r);
  449.  
  450. [K,mu]=h_inf(Gss,r,0,100,20);
  451.  
  452. ////////////SERE
  453. s=poly(0,'s')
  454. G=-4.714841d-3*(s**2+0.0218517*s+1.793965)*(s**2-0.3057194*s+4.752967);
  455. g=g*(s**2+0.4002371*s+4.812033)*(s-4.911808);
  456. g=g/((s**2)*(s**2+0.0221061*s+1.833885)*(s**2+7.392301d-3*s+3.18453));
  457. g=g/(s**2+2.716316d-2*s+4.624171);
  458. w1=0.33333333*s**3+0.3494322*s**2+9.157714d-2*s+0.012;
  459. w1=w1/((s**3)*(0.01*s+1));
  460. w2=2*s/(s+100);
  461. P=[w1,-w1*g;
  462.    0 , w2;
  463.    1, -g];
  464. [sk,mu]=h_inf(P,[1,1],1,1,1);
  465.